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Abstract. We describe an algorithm that takes as input a complex sequence 
(u n ) given by a linear recurrence relation with polynomial coefficients along 
with initial values, and outputs a simple explicit upper bound (v n ) such that 
\u„\ < v„ for all n. Generically, the bound is tight, in the sense that its 
asymptotic behaviour matches that of u„. We discuss applications to the 
evaluation of power series with guaranteed precision. 



1. Introduction 

A sequence u € C N is polynomially recursive, or P-recursive (over ( 
a non-trivial linear recurrence relation 



if it satisfies 



(1) pW (n) u n+s + • • • + (n) u„+i + p [0] (n) u„ 

with polynomial coefficients pW S Q[ti]- Likewise, an analytic function (or a formal 
power series) u is differentially finite, or D-finite, if it is solution to a non-trivial 
linear differential equation 



= 



(2) pW (*) «M (z) + • • • + pi 1 ! (z) u'(z) + pM (z) u(z) = 0, 



The coefficients of a D-finite power series form a P-recursive sequence, and con- 
versely, the generating series of a P-recursive sequence is D-finite. Numerous se- 
quences arising in combinatorics are P-recursive, while many elementary and special 
functions are D-finite. 



Starting with the works of Stanley (19801, Lipshitz (19891 and Zeilberger (19901 



D-finiteness relations have gradually been recognized as good data structures for 
symbolic computation with these analytic objects. This means that many opera- 
tions of interest may be performed on the implicit representation of sequences and 
functions provided by an equation such as ([T|, ( 2j) along with sufficiently many 
initial values (see Salvy and Zimmermann 1994| Stanley 19991. In recent years, 



significant research efforts have been aimed at developing and improving algorithms 
operating on this data structure. 

In this article, we describe an algorithm for computing upper bounds on P- 
recursive sequences of complex numbers. Specifically, we prove the following theo- 
rem (whose vocabulary is made more precise in the sequel). 

Theorem 1. Given as input a reversible recurrence relation of the form (fil with 
rational coefficients along with initial values defining a sequence (u n ) € , Al- 
gorithm^computes A G R+, K G Q, a G Q+ (the set of positive algebraic numbers) 



Key words and phrases. Algorithm, bounds, Cauchy-Kovalevskaya majorant, certified evalua- 
tion, holonomic functions. 
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and d> such that 



(3) VneN, \u n \ < An\ K a n (f>(n); 

with 4>{n) = e°^ n ' . Moreover, for generic initial values, k and a are tight. 

Asymptotic expansions of P-recursive sequences are a well-studied subject (see, 
e.g., Odlyzko 1995 Flajolet and Sedgewick 2009 1 and their computation has been 



largely automated (Wimp and Zeilberger 1985 Tournier 1987 Flajolet et al. 1991 



Zeilberger 20081. While an asymptotic estimate gives a precise indication on the 



behaviour of the sequence for large values of its index, it cannot in general be used 
to get an estimate for a specific value. Our result lets one obtain explicit bounds 
valid for any term, while the tightness of the bound with respect to the asymptotic 
behaviour implies that the bound is not straying too far away from the actual value. 
These bounds may be useful both inside rigorous numerical algorithms for problems 
such as D-finite function evaluation or numerical integration, or as "standalone" 
results to be reported to the user of a computer algebra system. The problem of 
accuracy control in several settings covering the evaluation of D-finite functions 



has been considered by many authors (see in particular Hoefkens 2001 Makino 



and Berz , 2003 Neher 2003 Rihm 1994 van der Hoeven 2003 2007). We review 



previous work on this problem in some more detail in S 5.2 Our main contribution 
from this viewpoint is to give bounds that are asymptotically tight. 

Example. To get a sense of the kind of bounds we can compute, consider the 
following examples. For readability, the constants appearing in the polynomial 
parts of the bounds are replaced by low-precision approximations, 
(a) Suppose we want to bound 



In = 



t n e 



2)I„ 



as a function of n e N. From the recurrence relation 2I n+3 = (n - 
and the initial conditions Iq,Ii,I2 < 1/5, Algorithm [5] finds that 

In < «! 1/2 2-" /2 • (0.26 n + 0.76) (" ^ 

In fact, /„ ~ n! 1 / 2 2 _n / 2_3 / 4 (7r/n) 3 / 4 as n — > oo, so that with the notations 
of Theorem 

« = 1/2, a = 2 1 / 2 are indeed recovered by our algorithm. 



(This example and the following one are adapted from Wimp and Zeilberger 
( |1985| Examples 2.1 and 2.3), who illustrate the computation of asymptotic 
expansions by the Birkhoff-Trjitzinsky method.) 
(b) The number t n of involutions of {1, . . . , n} satisfies the recurrence relation 

t(n + 2) = (n + l)t(n) + t(n + 1), t(0) = t(l) 
-1/4^,1/2^-1/4^-1/4 ag n 



1, 



Knuth 



1997 



§5.1.4). 



and t n - (87r)- 1 / 4 n! 1 / 2 eV"- 1 /4 n -i/4 as n ^ ^ ( see 
Assume that we wish to bound the probability that a permutation chosen uni- 
formly at random is an involution: the same algorithm leads tcQ 

^ < (0.90 n + 2.69) n^ 1 ' 2 [z n ] exp 
n\ 1 — 



= 0('< 



, 1 /4„|-l/2 p 2^¥i 



We use [z n ]f to denote the coefficient of z" in the power series /, see the end of §1 for 



notations. 
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Figure 1. Outline of our bound computation method. Solid ar- 
rows represent computation steps; dashed arrows indicate proof 
steps without counterpart in the algorithm. 



Compare (Flajolet and Sedgewick 2009 Example VIII. 5). Notice that, in 



addition to the parameters a and k of Theorem [T] the subexponential growth 
type e°^^i is preserved. However, our algorithm is not designed to preserve 
the constant in this O(-) term, 
(c) One of the fastest ways to compute high-precision approximations of tt resorts 
to the following formula due to Chudnovsky and Chudnovsky| ( 1988) p. 389): 

640320 3 / 2 , (-l) fe (6fc)!(13591409 + 545140134fc) 



DC 

E 

k=0 



12tt 



where ti 



(3/c)!(fc!) 3 640320 3fc 



Using the method of j |4.2| on the obvious first order recurrence relation satisfied 
by (tk), our algorithm leads to 



E 

k—n 



< 10 6 (2.3 n 3 + 13.6 n 2 + 25 n + 13.6)a™ 



where a 



0.66 ■ 10 



-14 



We see that each term of the series 



151931373056000 

gives about 14 more correct decimal digits of tt, and we can easily deduce a 
suitable truncation order to compute tt to any given precision, 
(d) Similarly, from the differential equation 

z Si"'(z) + 2 Si" (2) + z Si'(z) = 0, Si(0) = 0, Si'(0) = 1 

the result of our algorithm shows that the Sine integral special function may 
be approximated with absolute error less than 10" 100 on the disk \z\ < 1 by 
truncating its Taylor series at the origin to the order 74. 

Outline. Our approach is summarized in Figurejl] Consider a solution (u n ) of Equa- 
tion ([lj. Classical methods involving Newton polygons and characteristic equations 
allow to extract from the recurrence relation some information on the asymptotic 
behaviours that (w„) may assume. We use these methods to "factor out" the main 
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asymptotic behaviour, thus reducing the computation of a bound on \u n \ to that 
of a bound on a sequence of subexponential growth. This sequence is solution to a 
"normalized recurrence" computed in that step. Using the correspondence between 
P-recursive sequences and D-finite functions, we encode this sequence by a differen- 
tial equation satisfied by its generating function (^2j. Then we adapt the method 
of Cauchy-Kovalevskaya majorant series to bound this generating function. The 
key point here, in view of the requirement of asymptotic tightness, is to find a ma- 
jorant whose disk of convergence extends to the nearest singularity of the equation, 
thus avoiding the loss of an exponential factor usually associated with the majorant 
series method (Q. We show how to deduce several kinds of explicit bounds on u n 
and ^2 n u n z n from the asymptotic behaviour and the majorant series (Q- Finally, 
we introduce our implementation of the algorithms of this article and we briefly 
discuss their use in the context of high-precision numerical evaluation (^5|. 
Terminology and Notations. We let Q[n](5) be the algebra of recurrence operators 
with polynomial coefficients, viewed as noncommutative polynomials over Q[n] in 
the shift operator S : C N — > C N , (w n )n£N l— ¥ ( u n+i)nefi- Note that the sequences we 
consider are indexed by the nonnegative integers. Similarly, d stands for the dif- 
ferentiation of formal power series, and Q[z](9) for the algebra of linear differential 
operators with polynomial coefficients, written with d on the right. Noncommuta- 
tive monomials are written and represented in memory with the coefficient on the 
left and the power of the main variable S or d on the right. 

For any formal power series u £ C[[z]], we denote by u n (or sometimes by [z n ]u) 



the coefficient of z n in u. Following van der Hoeven ( 2003 1 , we also write 



DC 



^u k z k , u n . = u k z k . 

k—n k—0 



To avoid ambiguity, most other indexed names are writt en using bracke t ed su per- 



scripts, like p[°l in Equation ([lj. We use the notations of Graham et al. (19891 for 
the rising and falling factorials, namely x n — YYk=o( x + ^) an( ^ x ~ = Ofc=o — 

In the statement of algorithms, we employ expressions such as "set x > v'\ to 
mean "compute an approximation of v by excess (without any precise accuracy 
requirement) and assign it to x". 



2. Factorial and Exponential behaviour 

In this section, we collect classical results on the asymptotics of P-recursive 
sequences. These will both allow us to make precise statements about the tightness 
of the bounds we compute and serve as a guide to organise the computation in 
order to meet these requirements. Moreover, we state effective versions of some 
parts of the results, that constitute the first steps of our algorithm. 

2.1. The Perron-Kreuser theorem. A linear recurrence relation 

(4) pM (n)u n+s + ■■■+ pW (n)u n+1 + (n)u n = 0, 

or the corresponding operator 'J2p^S k , is called nonsingular when pW(n) ^ for 
all n £ N. It is called reversible when pt°' (n) ^ for all n £ N. 

Assume that the coefficients p[ fe l(n), k = 0, . . . , s of ([3} are sequences such that 
pM(ri) ^n^oo c k n dk for some cj~ £ C,dk £ Z (for instance, they are rational 
functions of n). If (u n ) is a solution of ([4} with u n+ i/u n ^ n ^oo then for 
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Figure 2. Newton polygons of recurrence operators, before and 
after normalization. 



the recurrence equation to hold asymptotically, the maximum value of dk + kn for 
k = 0, . . . , s must be reached at least twice, so that the corresponding terms can 
cancel. This means that — k must be among the slopes of the edges of the Newton 
polygon of the equation. 

The Newton polygon of Q is the upper convex hull of the points (k, dj.) £ K 2 , 
k = 0, . . . , s (see Figure |2j. If e is an edge of the polygon, we denote by — n(e) its 
slope. If (t, d t ) is the leftmost point of e, then the algebraic equation 

(5) Xe(A)= c kX k - l ^0 

is called the characteristic equation of e. Observe that the degrees of the charac- 
teristic equations sum up to the order s of the recurrence. 

Theorem 2 (Poincare, Perron, Kreuser). For each edge e of the Newton polygon 
of let Ae,i, A e ,2, • • • be the solutions of the characteristic equation \e, counted 
with multiplicities. 

(a) If for each e, the moduli |A e ,i| , |A e 2 1 ) • • ■ are pairwise distinct, then any solution 
(u n ) that is not ultimately satisfies u n+ i/u n ~„_>oo \ e ^n K ^ for some e and 
i. 

(b) If moreover Q is reversible, then it admits a basis of solutions (ut e,l ]) e! i<j< c i C g Xe 
such that 

[e,i] 

(6) Un+1 - A -n K(e) 

(c) If there exists e and i ^ j such that \\ e .i\ = |A e ,j|, results analogous to |ap and 
|ftp hold with the weaker conclusion 

(7) limsup 

n — >oo 

Definition 1 (Normalized Recurrences). If all the edges have nonnegative slope 
(i.e., if after dividing Q by pM, each coefficient tends to a finite limit as n — > 00), 
the recurrence is said to be of Poincare type. In that case, we call it (and the 
corresponding operator) normalized if the Newton polygon has a horizontal edge. 

Thus a normalized recurrence is one whose "fastest growing" solution has purely 
exponential (as opposed to factorial) growth. 



? !re(e) 



|A e 
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Algorithm 1: Factorial and exponential behaviour 



1 function Asympt(J2t=o blk] ( n ) Sk e QM^)) 

-1 dcg& w -dcgh l31 



s — k 

p a «- El= b d+tl zi where d = de § bls 

return (k, P a ) 



Item ([a]) above is known as Poincare's theorem (Poincare 18851; Items (|bj and 
(|c]| are Perron's theorem (Perron 19Q9a|b 1921 1 in the case of recurrence relations 
of Poincare type, and the Perron-Kreuser theorem (Perron 1910 Kreuser| |1914| in 
the general case. In addition to the original works, we refer to Meschkowski (1959) 



and Guelfond ( 1963 1 for accessible proofs of Poincare's and Perron's theorems. 



Various further extensions and refinements of these results are available, see, e.g., 



Schafke (19651, Kooman and Tijdeman (19901, Pituk (19971, Buslaev and Buslaeva 



(20051, and the references therein. 

In other words, the Perron-Kreuser theorem states that Q admits a basis of 
solutions of the form given by Theorem [2] in some neighborhood of infinity. The 
assumption that ([1} is reversible ensures that any solution near infinity extends to 
a solution defined on the whole set of nonnegative integers. 

2.2. Dominant Singularities. If P is a polynomial, we denote by ord(£,P) the 
multiplicity of £ as a root of P. We call dominant roots of P those of highest 
multiplicity among its nonzero roots of smallest modulus. We denote by <5(P) and 
ord^P) their modulus and multiplicity, respectively. By convention, the dominant 
root of a monomial is oo. We call dominant poles of a rational function the dominant 
roots of its denominator; and dominant singularities of a differential operator with 
polynomial coefficients the dominant roots of its leading coefficient. 

Besides standard symbolic manipulation routines, we assume that we have at our 
disposal a few operations on real algebraic numbers represented using the notation 
S(P), namely a function that decides, given F,Q G Q[z], whether S(P) < S(Q), 
S(P) = S(Q) or S(P) > 5(Q) and a procedure to compute arbitrarily good lower 
approximations of 6(P). The comparison can be based on a symbolic- numeric ap- 
proach as in (Gourdon and Salvyl 19961. Modern polynomial root finders such as 



MP Solve (Bini and Fiorentino 20001 or those of major computer algebra systems 
provide the required numerical evaluation features — and much more. Since we are 
interested only in S(P) as opposed to all roots of P, we may also use a simple 
procedure based on Graeffe's method (see, e.g. 



Schonhage 1982| §14) if no general 



polynomial solver is available. More generally, most steps of Algorithms [3] and [4] 
involving no precise accuracy requirement may be implemented using interval arith- 
metic or floating-point arithmetic with careful rounding instead of symbolically. 

Remark. Although we work over Q all along this paper for clarity, we expect that 
most results adapt without difficulty to any "sufficiently effective" subfield of C. 
However, the way to perform the basic operations we assume available in this 
section (as well as the details of some algorithms, especially Algorithm [3] below) 
may differ. 

2.3. Generic Growth of the Solutions. Let R £ Q[n](S) be a nonsingular 
reversible operator of order s. Then any solution of the recurrence relation R- u = 
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is uniquely determined by its initial values (uq, . . . , u s _i) G C s . Accordingly, we 
say that an assertion is true for a generic solution of R - u = 0, or for generic initial 
values, if it is true for any solution u such that (uq, . . . , u s _i) eC s \y where V is 
a proper linear subspace of C s . 

Theorem [2] implies that the factorial and exponential asymptotic behaviour of 
the "fastest growing" solutions is determined by the dominant singularities of R. 
We use Algorithm [l] to extract this asymptotic behaviour, which is in fact that of 
a generic solution of R ■ u = 0, as stated by Proposition [3] below. 

Proposition 3 (Factorial and Exponential Growth). Write R as X)fe=o b^ k '(n)S k G 
<Q[n\(S) and assume frMfrW for some k G {0, . . . , s — 1}. Algorithm^ computes 
(ft, P a ) — Asympt(TZ) such that for any solution (u n ) of R ■ u = 0, 

U n 1 /" , 1 

< a where a ■ 



lim sup 



the generic case. 



(8) 

with equality 

Proof. The inequality follows from Theorem [2] since —k is the slope of the rightmost 
edge e of the Newton polygon of R and P a is the reciprocal polynomial of \e- It 
remains to show that equality holds for generic initial values. Let V = ker R C C N . 
Also by Theorem [i] there exists G V such that 



lim sup 



1/r. 



a. 



This can be extended to a basis . . . of V. Let u = J2k ^' fe W fc ] G V. By 
construction of k and a, we have the inequality limsup \u n /n\ K \ 1 ^ n < a. Up to 
extraction of a subsequence we can assume (i) that Un does not vanish, (ii) that 



|iin'/n! K | 1 /™ — > a and (iii) that there exists /3 < a such that |w„/n! K | 1 ^ n 
n — > oo. Then 



/3 as 



aw + api^t + .-. + aW- 



1/n 



a 



so that /3 = a unless 



AP] U L 21 



-aw, 



which does not happen for generic AM. 



□ 



Accordingly tighter results hold if the assumptions of Theorem |2J|bJ are fulfilled. 

2.4. Generating Function and Associated Differential Equation. Consider 
again a nonsingular recurrence operator R = Ylt=0 °^S k € Q[n](5) (with b^°\ ^ 
0). Using the Euler derivative 9 — z4-, it is classical that the generating series u{z) 
of u G keri? cancels the associated differential operator D = Ylk=o a ^9 k G Q[z](0) 
computed by RecToDiffeq (Algorithm plrl Dividing out by o", this rewrites 



(9) 



1 



[o] 



0. 



2 Actually, the classical translation of recurrence operators to differential operators uses 9 = 1. 
The multiplication by g in our version comes from our choice to use sequences indexed by N rather 
than Z. 
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Algorithm 2: Recurrence to normalized differential equation 
1 function RecToDiffeq(i? = Efc=0 &Ws * G QN(S')) 



g <- n/gcd(6W,n) where n = IlLiO 

expand Eft=oEj c fe^ 
return Z) 



s-ki 



- k) j asJ} = Efc=o oI * le * 



k) 

>thus R=YZ=oCkjS k (n-ky 



6 function Normalize(i? £ Q[n](S), k £ 



pi q < — k (in irreducible form, with (p, q) = (0, 1) if k = 0) 
compute the symmetric product R — EfcLo & k \n)S k of i? and 
(n + 

> see, e.g., 



1 



StanJey (1999 §6.4) 



return RecToDiffeq(i?) 



A point zq £ C is a regular point of ([9]) if any solution it has polynomial growth 
u{z) — l/\z — Zq\° as z — » 2o m a sector with vertex at zq. Regular points 
encompass ordinary points, where the equation is nonsingular and thus has analytic 
solutions by Cauchy's theorem, and regular singular points. Fuchs' criterion (see, 



e.g., Ince 1956 §15.3) states that is a regular point if and only if for all k, the 
coefficient a^/a^ of ([9} is analytic at 0, while zq ^ is a regular point if and only 
if each a^/a^ has a pole of order at most r — k in Zq, (This criterion still holds if 
the a^l/aM are replaced by meromorphic functions.) 

Lemma 4. If R is normalized (Definition^, then the origin is a regular point of 
D, and the reciprocal polynomial of the leading term w- r ' of D is the characteristic 
equation of the horizontal edge of the Newton polygon of R. 

Proof. Using the notations of the function RecToDiffeq() in Algorithm]^ let = 
degfrM f or a n anc j m = degg. Thus r — max|_ S k ' + m. The leading term 
of z~ k as an operator in 9 with Laurent polynomial coefficients is z~ k 9 : > , hence 
a\ r \z) = ~J2k=o c kr zS ~ k - The condition that R is normalized translates into d^ = 
maxj^Q rfW, that is, = d^ = r — m for some k < s. It follows that a^(0) = 
c sr ^ 0, hence is a regular point by Fuchs' criterion. Finally, if R is normalized 
and if e is the edge of its Newton polygon such that n(e) — 0, then the general 
expression 

Xe(A) = A~* a k,dW^ k 

d lk] +kK(e) 
=d [s '+SK(e) 

(where t is such that x e (0) ^ 0) simplifies to % e (A) = A~* EdM =r a k,r^ k - O 

In the general case, we normalize R by a change of unknown sequence preserving 
P-recursiveness before we compute the associated differential equation. This is 
described in the next proposition. Figure [2] gives an example of normalization of 
recurrence operators and of its action on their Newton polygons. 

Proposition 5. Let R £ Q[n](S) be nonsingular, reversible, with nonzero con- 
stant coefficient with respect to S. Let (p/q,P a ) — Asympt(i?) as computed by 
Algorithm^ and assume that S(P a ) < oo. Algorithm^ computes a normalized 
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differential operator D — Normalize(i?,p/q) that cancels u(z) — X)^Lo ' l PnUnZ n for 
all sequences ip and u such that 

(n + q) p ip n + q = ipn and R ■ u = 0. 

The origin is a regular point of D, and the modulus of the dominant singularities 
of D equals d(P a ). 

Proof. Let a — l/5(P a ). Let (u^\ . . . , uW) be a basis of keri? having the asymp- 
totic behaviours given by ([7f . In particular limsup n _ >00 |u'' ! ] /n\ p / q \ 1 / n < a for all k. 
Let , . . . , ■01' 3-1 ]) be the basis of solutions to (n+q) p -i/Vi+g = V'n corresponding to 
the initial values ip^ = Sij for < i, j < q, where 5ij is the Kronecker symbol. Algo- 
rithm [2] constructs R such that for N large enough, the sq sequences {4>n^u$) n >N 
generate {it \ (R ■ u) n = for n > N}. For all j and k, limsup|'0n'wn'| 1 ^ n — a - 
Assume that u — k ^'^ip^u^ is solution to R ■ u = in some neighborhood of 
infinity. Then limsuplitnl 1 /" < a (indeed, if e > 0, then |m„| < (a + e) n for n large 
enough). On the other hand limsup|'0Ti''Un' I 1 /™ = a for at least one {j,k). Hence, 
by Theorem [2] the operator R is normalized and the largest modulus of a root of 
the characteristic equation associated to the horizontal edge of its Newton polygon 
is a. Applying Lemma [4] concludes the proof. □ 

In the sequel, we will choose as normalizing sequence the solution to (n + 

q) p ipn+ q = i> n g iven b Y 

^ n = q->T{n/q + l)-v. 

Observe that (VvOneN is monotone: indeed, the function x 1— > g x r(a;-|-l) is increasing 
for x > as soon as logg > 7 (the Euler-Mascheroni constant), and the remaining 
case q = 1 is obvious. 



3. SUBEXPONENTIAL BEHAVIOUR: MAJORANT SERIES COMPUTATION 

The results of the previous section allow us to compute the generic factorial 
and exponential asymptotic behaviour of solutions of a linear recurrence relation 
with polynomial coefficients. We now turn to the computation of a bound for the 
remaining subexponential factor of a particular solution. 

3.1. Majorant Series and the Cauchy-Kovalevskaya Method. The main tool 
we use is a variant of the Cauchy-Kovalevskaya majorant series method, which 
usually serves to establish the convergence of formal series solutions to differential 
and partial differential equations, but may also be applied to obtain explicit bounds 



on the tails of these solutions (see also 6 5.2 for more on this 



Definition 2 (Majorant series). A formal power series v S R+[[z]] is a majorant 
series of u G C[[z]], and we write u < v, if \u n \ < v n for all n£N. 

In particular, the disk of convergence of v is contained in that of u, and if z lies 
inside the disk of convergence of v, we have that |w n; (z)| < i>„ ; (|z|) for all n > 0. 
Other immediate properties of majorant series are summarized in the following 
lemma. 
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Algorithm 3: Tight majorant series for rational functions 

function BoundRatpoly(r = N/D € Q(z),P a G Q[z],m€ N*) 
let A + B/D = r with A,Be Q[z}. 

compute the squarefree factorization D = D\D\ ■ ■ ■ Di of D 
compute the coefficients h^d € Q[C] of the partial fraction decomposition 

B(Z) _ ^pfc yii fej,d(C) . C pp p p. 

D(z) ~ Z-f»=l ^15,(0=0 l^d=l ((-z) d " oee > 

for i = 1 , . . . , fc do 
for d = 1, . . . , i do 

MOM 

set Nq > max(l, 1 + deg A, max*; 



Bronstein 



(2005 



§2-7) 



=m+l log(5(D,)/<5(P»))/ 



d=0 



1 <^^£ (S(P a )/S(D t )Y 



'(n+iy 



No-1 



lett(n)=E-=iE 

compute the truncated series r.^ g (z) = Y^rl= 
set h(N ) >max^ 1 (|r„|/(("+™7 1 )J(F Q )«)) 
return an approximation by excess of max (/i(7V ), i(iVg)) 



Lemma 6. Assume that u, it^ , vffl <E u,i>M,?;[ 2 l 6 are stte/i 

it < u, uM < and iJ 2 ] < I;! 2 !. Then 

tt [il + „Pl< u H] +fl P!. tt [i] u [a] < v m v m. tt M ou [iI< e P] OB [i] 

dz dz 

where in the last inequality it is assumed that uM(0) = i>M(0) = 0. 

In the neighborhood of an ordinary point, majorant series for the coefficients 
of a differential equation like Q give rise to similar majorants for the solutions. 
Indeed, if 

{ „M = 6l-V^) + • •• + ^ |U(0)I " U(0) ' ' (0)l " V (0) 

where a[ fe ],&[ fc J are analytic functions at such that < for all fc, then by 
induction u < u. This result does not hold if one of the a'*' has a pole at 0; however, 
the method may be adapted to the case where is a regular singular point of the 
differential equation. We give one way to do this in §3.3| for a more complete intro- 
duction to the "usual" Cauchy-Kovalevskaya method in the ODE setting covering 



the regular singular case, see Mezzino and Pinsky 



statement along these lines, see van der Hoeven (2003 Proposition 3.7). In any 



(19981, and for a more general 



case, the first step for obtaining majorant series for the solutions of a differential 
equation using the Cauchy-Kovalevskaya method is to compute majorants for its 
coefficients, which in the case we are interested in are rational functions. 

3.2. Bounds for Rational Functions. Consider a rational function r(z) = N(z) / D(z) 

r nZ n , D(0) 7^ 0. The sequence (r„) satisfies a linear recurrence relation with con- 
stant coefficients, whose characteristic polynomial is the reciprocal polynomial of 
D. This recurrence can be solved by partial fraction decomposition of r, yielding 
the explicit expression (recall that x- and x n denote respectively the falling and 



EFFECTIVE BOUNDS FOR P-RECURSIVE SEQUENCES 



11 



rising factorials) 

ord(C-D) 

(10) r n = £ /i [c ' <fl -(n+l)* zr -C~ n , n>max(0,deg7V-degD+l), 

D(()=0 d=l 

with h^' d ^ S Q(C)- We are now aiming at a bound of the form |r„| < M5(D)~ n n mds D . 
In view of later needs, Algorithm [3] takes as input a polynomial P a and a positive 
integer to. It returns a bound of the form r(z) < M(l — az)~ m , where a = l/S(P a ). 
In particular, when P a — D and to = ords(D) this bound is tight. 



To compute a suitable M, we start with the right-hand side of ( 10 1 divided by 

b n = [z n ] jz — ^r- = (n + l)"^ 1 • a n . 
1 l (l-az) m y ' 

By applying the triangle inequality, we get a sum t(n) of terms of the form 

(n+ l) m_1 

where < c, < A < 1, and m < d only if A < 1. Such a term is decreasing for 
n > 1 if <i < to and for n > (d — to) / log(l/A) otherwise. We compute an index iVo 
starting from which the inequality |r n /6„| < t(n) is guaranteed to hold and t(n) 
is guaranteed to be decreasing; then we adjust M from the explicit values of the 
first Nq coefficients and bounds on the tails. 

For this last part, consider the squarefree decomposition D — D\D\ ■ ■ ■ D\. If £ 
is a root of Di, then each h^' d ^ may in fact be written h^ ,d ^ = hi t d{C) ■ C,~ d for some 
polynomial hi 4 E Q[(] depending only on Di and d. Moreover, in this expression, 
ICI 1 may be bounded by 5(Di)~ 1 . Hence we have 



(11) 



*~ n j2 J2 E^(0c d(n+l) "~' ' 

i=l Di(()=0d=0 





We may take for t(n) the right-hand side of ([TTJ, or even a suitable numerical 
approximation. To deal with the sum in parentheses, we may bound C _d ^i,d(C) 
term-by-term, replacing once again ( e by 8(Di) 1 or 6(( degDi Pi(l/Q)~ e depending 
on the sign of I. We may also simply compute low-precision enclosures of the roots 
of Di and then use interval arithmetic. 

The complete procedure is summarized in Algorithm [3] We have thus proved 
the following. 

Proposition 7. Given r = N/D € Q(z) (in irreducible form) , P a G Q[z], and to € 
W, such that < 8(P a ) < 5(D) and 8(P a ) = 5(D) only ifm> ord s D, Algorithm^ 
computes M = BoundRatpoly(r, P a , to) £ Q + satisfying r(z) < M(l — z/S(P a ))~ m . 

To improve M, we may loop over lines [10] and [TT] of Algorithm [3] doubling N 
each time, until Nq or t(N ) — h(N ) reaches some specified value. 
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Algorithm 4: Majorant series for normalized D-finite functions 



1 function BoundNormalDiffeq(X^ = o °> [k] k G Q[z]{0),P a £ 



\z\,u..) 



for k = 0, . . . , r — 1 do 

cl fc l <— (a[ fe '/ a ^)z=o ( or fail with error "0 should be a regular point") 



( a [ fc ]/a [rl 



T <- max{0;ord 5 (denaL fe J) -r + fc | < k < r - 1 and 8(dena^) = S(P a )} 
for k = 0, . . . , r — 1 do 

MW <- BoundRatpoly(aW , P a) T + r - k) > thus 
aW <M[ fe ](l-az)- T " r+fc 

M^max^MW/C"" 1 ) 
compute JfeN' such that K > 2M5(P a ) 
starting with N2 = 1, double iVa until 
Y,l~=lW k] \ N 2<(l-M6(P a )/K)N% 

compute u-n 2+ i and v-n 2 +i where v is given by (18 1 with A = 1 
A <- max^o |tt n | /v n 



return (T, if, A) 



3.3. Bounds for D-finite Functions. We now apply the Cauchy-Kovalevskaya 
method to deduce a majorant series for u(z) from the asymptotic behaviour of (u n ) 
obtained in §1 and majorant series for the coefficients of an associated differential 
equation. The majorant series we obtain is "simpler" than u(z) in the sense that it 
always satisfies a differential equation of order 1. 

By Fuchs' criterion, we may isolate the constant term of each coefficient of ([9}, 
giving 



(12) 



Q(9) 



z(a 



x\r—l\ nr—1 



+ ar°]) 



where Q £ Q[X] is a monic polynomial of degree r and the aM are rational functions 
of z. Let rrik £ N be the maximum multiplicity of a point of the circle \z\ = S(P a ) 
as a pole of 5,^ and let T = max(0,max^J(mj. — r + k)). We emphasize that, 
although Algorithm |4]takes P a as input, the whole point of the method is that S(P a ) 
may indeed equal the modulus of the dominant singularities of D. In that case, the 
integer T is sometimes called the Malgrange irregularity of these singularities (see 
Malgrange 1974 1, and by Fuchs' criterion again, T = if and only if the dominant 



singularities are all regular. Using Algorithm [3] we compute bounds of the form 



(13) 



aw < 



(1 - azf- k + T 



al fc i 



< M [k] 



-k+T-1 
k + T-1 



for the coefficients of the equation, with a = l/S(P a ) as usual (lines |6fj7] of Algo- 
rithm |i} . 



Extracting the coefficient of z n in (12 I, we get 



(14) 



n— 1 r— 1 

Q(n)u n = J2Y,~ a n-i-ji ku i- 

j=0 k=0 
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Since Q is monic, let N% be such that Q{n) > for n > Ni; then by (13), for such 

n, 



(15) Q(n) \u n \ <Y.Y. M[k] 

j=0 k=0 



n—l—j + r—k+T — 1 
r-k+T - 1 



„ n—l—j -k I I 

a J j \uj\ 



Lemma 8 (Reduction from order r to order 1). Let M = max^J /( r fc 1 ) and 
< j < n — 1; then 

g A/W ^-l-J + r-k+T-l^ jk ^ Mn r-i fn-l-J+T 



fe=0 



r-k+T-l 



T 



Proof. For k < r — 1, we have 

n—l—j+T\ 1 / n—l—j + r— k+T— 1\ (n — j+T) 



T 



r-k+T-l 



\r-l-k 



(T + l) 



r — 1 — fc 



< (n - j) 



r-l-fc. 



thus 



rc-l-j+T' 
T 



-1 r-l 



k=0 



-l-j + r-k+T-l\ k 
r-k+T-l J 



<£ M M/ (n _ jT -i-< 
< Mn r -\ 



fe=0 



establishing the lemma. 



□ 



With M as in Lemma [sj choose K > M /a. Let N 2 > Ni be such that Mn r < 
aKQ(n) for n > N%. Suppose that some sequence (v n ) satisfies v n > \u n \ for 
< n < N 2 and 



(16) 



1 71 — 1 

v n = -T i K 



3=0 



n-l-j + T 
T 



for all n > I. Let n> N 2 . Assuming \v,j | < Vj for all j < n — 1, and using ( 15 I and 
Lemma [8] we get 



T 



hence by induction |u„| < w„ for all n G N. Now (16) translates into 



(17) 



(1 - m) 



T+ 



T«( 2 ): 



which admits the simple solutions ( 18 1 below. 

Finally, we adjust the integration constant A so as to ensure that \u n \ < v n for 
n < N2 (lines 11 12). If no specific solution of ([9} is given (i.e., if we drop the 
parameter u. n of Algorithm [Ij we still obtain a result valid up to some multiplica- 
tive constant by simply ignoring this last part. The result of this computation is 
summarized in the following. 

Proposition 9. Let D S Q[z](9), and let u :n be a function that computes truncated 
series expansions of a specific u € kerD up to any order n. Let P a £E Q[z]. Assume 
that is a regular point of D and that the dominant singularities of D are finite 
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Algorithm 5: Bounds for general P-recursive sequences 

1 function BoundRec(i? = EjUo^W 5 * e Q[n](S), [it Q , . . . ,u 8 _i] G Q[i] s ) 
R- S- m where m = min{fc | ^ 0} 
(k,P«) <- Asympt(i?) 

> Normalize and encode the subexponential part by a differential equation 
D <— Normalize (i?, k) 

> Bound the solutions of the differential equation 
define a function u- : . that "unrolls" the recurrence relation R ■ u = 
starting from Uq, . . . , it s _i to compute 

";n = ELo q~ pk/q T(k/q + l)~ p u k z k (where p/g = k) for any n £ N 
(T, K, A) <- BoundNormalDiffeq(£>, P a ,u..) 
return (k, T, P q , if, A) 



and of modulus at least 5{P a ). Then BoundNormalDiffeq(Z), P a . u-.) (Algorithm^ 
returns T £ N, K £ W , A £ Q + such that 



(18) u(z) < v(z) = 



- (1 



K / T 

A exp -= otherwise. 

(1 — az) 1 



In addition to its modulus a, Algorithm 4 actually preserves the irregularity T 
of the dominant singularity of the differential equation, which is connected to the 
subexponential growth of the coefficient sequence. 

Remark. Sometimes all we need is a simple majorant series satisfying the tightness 
property of Theorem [I] for the solutions of a differential equation of the form ^ 
at an ordinary point. Instead of the results of this section, we may then apply the 
"plain" Cauchy-Kovalevskaya method outlined in |3.1| using a majorant equation of 
the form 



(l-az) N ^\ k J Vl 

v ; k=0 v 7 

This gives the majorant series v(z) = exp(M/(l — az) ). If additionally the 
dominant singularity is regular, we may instead use the Euler equation 

^ AfM 
(1 - az) r - k ' 



^ (i 

k=0 



yielding v(z) = A/(l-az) A where a r \ r - M^-^a^ 1 X^ 1 = 0. In both 

cases suitable parameters M, resp. may be determined using Algorithm [3] 

4. Explicit Bounds 

4.1. P-Recursive Sequences. At this point, we are able to bound u n by a se- 
quence v n given by its generating series v(z) — C' pq v(z), where v is an explicit 
series satisfying a differential equation of the first order, and we have denoted 

00 

<vw -££»■■ 

n=0 
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(Note that series whose coefficients satisfy recurrence relations of the first order, 
that is, hypergeometric series, cannot serve as asymptotically tight bounds for 
normalized D-finite functions because the range of asymptotic behaviours that their 
coefficient sequences assume is not wide enough: their "subexponential" asymptotic 
growth is always polynomial.) 

Proposition 10. Given as input a nonsingular reversible recurrence operator R G 



[n] (S) along with initial values uq, 



-l G Q[i] defining a solution (u n ) G 



of R ■ u — 0, the function BoundRec (Algorithm^ computes p/q G Q, P a G 
T G N and K, A G R+ such that 



(19) 

where v r< 



Vn G N, |ttn| <«„ = g?"r(- + l) 



1 V„ 



defined as in (18 1. Additionally, for generic (uq,. 

I In 

am sup — 



= 1. 



Allowing initial conditions in Q[i] rather than Q is convenient in view of some 
applications to numerical computations with D-finite functions (Sp). 



Proof. This follows from combining the statements of Propositions [3] [5] and |j] 
Recall that we have chosen ijj n — q~v n T(n/q + l)~ p . After Linejijof Algorithm [5] 
the operator R satisfies the hypotheses of Proposition [5] Hence the operator D 
computed on Line [4] cancels u(z) — X)^Lo i } nU n z n , and the function u-. defined on 
the next line does indeed compute truncations of this series. By Proposition [9] it 
follows that u < v and, multiplying the coefficients by if}~ 1 , that u < v. Finally, for 
generic initial values, 



lim sup 

n — >oo 

by Proposition [3] 



l/r. 



lim sup 



,|K ,n+o(l) r) 0(l) 



l/r 



□ 



Although this representation (19 1 is satisfactory for many applications, more 
explicit expressions for the coefficients v n are sometimes desirable. If T = 0, it is 
readily seen that 



(20) 



= At 



n + K - 1 
K - 1 



For T > 0, the general coefficient v n still admits a rather complicated "closed-form" 



expression in terms of the general hypergeometric function F (see Graham et al. 
1989 §5.5): one may check that 



Aa' 



OO 

E 

k=0 



1 (Tk 



- n - 
n 



( n+T 

Aa n tFt \ t+i 



n+T+1 

T 
T+2 

T 



K+2T-1 
2T 



However, v n may in turn be bounded by much simpler expressions without losing 
the asymptotic tightness (in the se nse of Theorem [T| using a simple version of 



the saddle point method (see, e.g., Flajolet and Sedgewick 
\z\], for any t G (0;l/a), we have v n < v(t)/t 



2009 



§4.3). Since 

v G K + [|zJJ, for any t G (0;l/a), we have v n < v{t)/t n . For fixed n, the right- 
hand side is minimal for the unique t n G (0; 1) such that Kat n = n(l — at n ) T+1 . 
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Asymptotically, t n satisfies 1 — at n ~ (K/n) 



1/(T+1) 



as n — ► oo. This approximation 



suits our purposes well: indeed, we set 

1 

r n = - 
a 



(21) 



1 



K 



(The term K + 1 in the denominator does not change the asymptotic behaviour 
and is such that r n £ (0; 1/a).) For T > 0, we obtain (with A = 1) 



w(y n ) 



(22) 

and similarly 
(23) 



a 



( K ) 
\n + K+\) 



exp 



K fn + K+ 1\ T+1 
T 



K 



. n /n + A+l\*"/ 



A' 



n + A+ 1 



a™expO(n 



T/(T+1)\ 



if T = 0. 

Going back to u„ itself, ( 22 I and ( 23 1 extend to bounds of the form ^ , that make 
the asymptotic behaviour u n — nl K a n apparent, by means of the following 
relation between ip n and n\ K . 

Lemma 11. For q £ N \ {0} and n > 3q/2, 



Tin! a + W« » < / {27r)P/q {n/q + 1)P P > ° 

'■ ~~ l T V9 + L ) Q ^ S —p/q„W/q - 



p < 0. 



Proof. Since r(x) is increasing for x > 3/2, 



9-1 9-1 

J] r(n/g + k/q) < T(n/q + l) q < ]J T(n/q + k/q + 1). 
k=0 k=0 



By Gaufi' multiplication theorem (see Abramowitz and Stegun 1972 Formula 6.1.20) 



9-1 



T(qz) = (2n) ( - 1 -^/ 2 q^ 2 l[T 



z + 



(Z £ C), 



fe=0 



this implies that 



( 27r )(9-l)/2 g"r(n/g+ 1)1 ( 27r )(9-l)/2( n+ !)9-l 

nq-V2 - T(n + 1) ~ g?- 1 / 2 

and the result follows by raising either inequality to the power oip/q depending on 
the sign of p. □ 

This concludes the proof of Theorem [l] 

Remark. If we content ourselves with computing a numerical bound for one coeffi- 
cient (or one tail, see next section) of a D-finite power series — that is, a bound for 
fixed n, as opposed to a formula giving a bound as a function of n — then majorant 
series with the same radius of convergence as the coefficients of the equation (and 
thus the method of 1 3.3 1 are not strictly necessary for the bound to become ulti- 
mately tight as n approaches infinity. Consider for instance Equation (JlJ in the case 
where is an ordinary point, and assume v > a with the notations of p.3| |Van der| 



c 
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Hoeven| ( |2003| §3.5) proves that if /p M < M(u)/(1 - vz) for k = 0, . . . , r - 1, 
then 

C 

U( - Z > ~ (1 - VZ )\{M(v)+l)/^ 

where C does not depend on v. Also assume that the majorizing procedure for ra- 
tional functions used to compute M(y) is tight enough to ensure that M(v) — 
0(n d (a/v) n ^ (as is Algorithm |3j with d = max^Jm/j). In a manner some- 
what reminiscent of the saddle-point method, we then choose, say, v — v n — 
(1 + l/n 1 /( 2d ))a, hence getting 



n-\-n 



\u n \ < v n = a 

This suggests that it is sensible to take v = (1 + l/n @ ^ 1 ' d ')a in the algorithms of 
van der Hoevenl (|2001| 12003 1. 



4.2. Tails of Power Series. In Examples |lj|c| and (|dj, the sequence for which 
we compute an upper bound is the tail t n = u n -(l) of a convergent series whose 
coefficients u n are given by a linear recurrence relation of the form (JTJ. In such 
a case, the sequence t n is also P-recursive, but its initial values are unknown — if 
we have in mind the evaluation of the sum of the series, these initial values are 
precisely what we are after. However, if u(z) < v(z), the general properties of 
majorant series (§^| ensure that |u„ ; (l)| < v n -(l). To avoid repeated majorant 
computations when working with D-finite power series, notably in the context of 
numerical analytic continuation (see |5.2| , we actually consider the slightly more 
general problem of bounding the tails u„ ; (z) of the j-th derivative of u at any point 
z such that \z\ < 5(p^), where i s the leading term of a differential equation 
with polynomial coefficients annihilating u{z). 

We assume once again that we have computed k = p/q and v such that u(z) < 
v(z) = C' p q v(z) (with p < 0, so that the radius of convergence of v is positive) 
using the algorithms of §S]and §S] The letters a, T, K denote the parameters of v 
appearing in ( [T8| . The formalism of majorant series proves handy here, as we have 
(z)| < UnP(|^|) by Lemma ^1 Notice that if p < 0, the point z lies within the 
disk of convergence of v but not necessarily in that of v, 

Proposition 12 (Bound on u n .(z) for large n). With z and v as above, assume 
that 



(l-alzD-'-'K, 



K=0 



(24) n > 



{a\z\)-*l* 
Then for all j, we have 



K 



(25) 



< 



(a\z\)-i/P + K + 1 
v U \r n ) M\ n ,(\z\ 



, K < 0. 



g -i"r(| + i) 



—( B )h(m, 

l)-p\r n J \r n J 



where r n is given by (21 1 and 

<?-i 



1 

= z 77 s — y^x u (= 1/(1 -x) for k = 0, i.e. p/q = 0/1). 



u=0 



The bound (25 I is generically tight up to subexponential factors. 
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Figure[3]illustrates the behaviour of this bound for entire functions, in the typical 
situation where the Taylor series at the origin "starts converging" only beyond a 
significant "hump". Once again, the factor n\vli in p5| can be brought out explicitly 
if desired using Lemma [TTj 



Proof. In the case k = 0, the condition (24 1 ensures that \z\ < r n < a^ 1 . Using 
the relation v n = ip n v n and the saddle-point bound % < v(r n )/r^ (notice the n), 
we obtain 

Yn x 1 n ' , n "Vn+k x 1 n ' 
k—O 

This proves ( |25| for k — 0. Now assume p < 0, and recall that in this case 
i/} n — q~ p / q T{n/ q + l) _p is increasing: hence 



x 



tq 



oo , cjo q— 1 , q — 1 oo 

fc =o ^ n+k t=o u=o «=o t=o ((n + <?) (n + 2g) . . . (n + t<?)) 

for n > x~ q / p . But this last condition follows from ((24ll since 



< h(x) 



\z\\-q/p 



<{a\z\)-*l> 



(a|z|)-«/p + J ftT+ 1 



<///' 



as soon as n > (a \z\) q / p , itself implied by (24 1. 

The estimates p2|), (23 1 still hold, hence the tightness of the bound. 



□ 



Bounds on u n .(z) are sometimes useful also when the condition (24 1 fails to be 
satisfied, especially for n — 0. Simple bounds independent on n give good results. 
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Proposition 13 (Bound on u n; (z) for small n). For all n £ N and < r < a 1 , 

r vW(\z\) k = 

(26) 



E(~f «< 



0. 



U=0 



Proof. The proof is similar to that of Proposition 12 For k = the result is 
obvious. Assuming n < 0, it holds for all x > that 



9-1 oo t «-l oo (_P x -q/p)-pt 

E^E*" E J:<E^ E q } 



k=n r " u=0 t=|n/?J Tqt «=0 t=L"/<?J 
since W = q- pt tl~P > [,-q/p)~ vt [.-pt)\ (t £ N); whence 



(-Pt)l 



In the important case where k = T = and K £ N, the u n; (z) actually admit 
closed-form expressions of the form (az) n p(n), wherep £ Q(az)[n]. Indeed, starting 

E^iC [il (rc)0 + l)^ we get 

K 



from (18 1 and writing (for fixed K) (n+fc + 1) 



if— l 



1 



(1 — QZ) 



A' 



(^-1) 



TvE( 



fe=0 



v ' (K-l)\^(l-az\ 



c W (n). 



This is the kind of formula that appears in Example flflc| . Such bounds are easier 



to read than (25 1, but they are numerically unstable due to cancellations. In a 



system providing numerical routines for hypergeometric functions, one can use the 
alternative expression 



(l-az) K 



{az) n 



K - 1 



QZ 



which does not suffer from this shortcoming. 

Finally, note that it might be worthwhile looking for refined bounds in applica- 



tions where T is large and \z\ ~ a -1 , since (25) becomes tight only for very large n 
in this case. Similar issues exist when K is too large; they may be mitigated by mod- 
ifying Algorithm [3] to compute bounds of the form p(z) + M/(l — az) m ,p £ Q+[z], 
which allows for a tighter choice of K. 

5. Applications and Experiments 

5.1. Implementation. We have implemented the algorithms described in this ar- 
ticle (with slight variations) in the computer algebra system Maple. Our imple- 
mentation is part of a submodule called NumGfun of the Maple package gfurj^] but 
the code computing bounds is largely self-contained. It provides routines that com- 
pute majorant series for rational polynomials (following |3.2| and D-finite functions 



(£ 3.3 £4.1 ) , and symbolic bounds for P-recursive sequences specified either using 
recurrence relations (S 4.1 1 or as tails of D-finite series ( |4.2) . All examples of this 
article were computed using this implementationF] 



"http : //algo . inria.fr/libraries/papers/gfunJitmlJ 

4 To be precise, using gi'un v. 3.48 under Maple 13. 
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It is also used by the Dynamic Dictionary of Mathematical Functions^ an inter- 
active web-based handbook of D-finite functions currently under development. All 
contents of the Dictionary are automatically generated from a compact description 
of each function (basically, a differential equation and initial values) using a mix 
of symbolic computation algorithms and document templates. The webpages the 
system produces are interactive in that they allow the user to trigger more com- 
putations, typically by asking for "more terms" in an asymptotic expansion. This 
is a situation where being able to display human-readable formulae rather than 
merely computing numerical bounds represents a significant benefit. Code based 
on this article provides majorant series for the Taylor expansions of the functions, 
truncation orders for these expansions to reach a given accuracy over a given disk, 
and symbolic bounds for their tails involving the truncation order. 

5.2. Application to the Numerical Evaluation of D-Finite Functions. Guar- 
anteed numerical computation with entire classes of functions usually involves the 
automatic computation of error bounds relating approximations, e.g., by truncated 
power series, to the functions they approximate. Elementary results from real and 
complex analysis commonly used to compute such error bounds include the alter- 
nating series criterion, Cauchy's integral formula, and several variants of Taylor's 
theorem. Karatsuba describes algorithms with error bounds for the evaluation of 



various special functions, including the hypergeometric function 2-P1 (see Karat 



suba 1999 and the references therein). Du and Yap (20051 provide bounds for the 



tails of the general hypergeometric series, where the parameters are allowed to vary, 
based on a detailed analysis of the variations of the coefficient sequence. For the 
more general case of D-finite functions, another ad hoc method is given by |van der| 
Hoeven (19991. In a different context, Neher (2003) uses Cauchy's estimate and 



complex interval arithmetic to bound the coefficients and tails of series expansions 
of arbitrary "explicit enough" analytic functions. This method is implemented in 



ACETAF (Eble and Neher 20031 



A further classical tool is the Cauchy-Kovalevskaya majorant series method dis- 



cussed in £3.1 This idea is exploited by van der Hoeven (2001 §2.4) to bound 



the tails of power series expansions of D-finite functions in the neighbourhood of 
an ordinary point of the equation, and later again in a much more general setting 



covering a wide range of functional equations ( van der Hoeven 2003 1 . This is the 



approach we rely on in this article: indeed, the algorithm we described in §!3.3| may 
actually be seen as a refinement of those suggested in §3.5 and §5.2 of the latter 
article. The main originality of our approach is the asymptotic tightness of the 
bounds. 

Finally, it should be noted that in the context of numerical evaluation, instead 
of using a priori bounds, it is often easier to compute successive error bounds in 
parallel to successive approximations of the result, until the desired accuracy is 
reached. The computation of validated numerical enclosures of solutions of ODE, 
DAE and more general functional equations has been the subject of extensive liter- 



ature since the sixties (see Rihm 1994 1 in the area of interval methods. Of special 
interest when working with power series is the integration of differential equations 



using Taylor models (see Hoefkens 2001 Neher et al. 20071. Taylor models are 
one among a fair number of different symbolic-numeric representations of functions 



http : //ddmf .msr- inria . inria.fr/ 
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used in interval arithmetic, several of which have a similar approach of bounds for 
solutions of functional equations: for more on Taylor models and their relation to 



other interval methods, see (Makino and Berz 2003 Neumaier 2003). Some of 



these methods were imported to computer algebra and revisited by van der Hoeven 



(20071 in the context of rigorous effective complex analysis. 



In a nutshell, the common idea is to write the (differential, say) equation at hand 
in fixed-point form u — <f?(it), where <I> is an integral operator, and to consider the 
action of $ on truncated power series augmented with error bounds, using rules 
such as 

j (a + ait + a 2 t 2 + [a,(3])dt C J ( fflQ + a\t)dt + -B^sy) + [a, 0\ • B{x). 

Here B(p) is an interval containing the range of p{x) obtained from the range of 
x. One then computes an approximate solution in the form of a Taylor expansion 
p(x) = a + ■ ■ ■ + a n x n and iteratively searches for a tight interval [a, /?] such 
that <I>(p + [a, f3]) C p + [a, (3), possibly narrowing the range of x or increasing the 
expansion order n as necessary. Under mild assumptions, the existence of such 
p + [a, (3] implies that of an actual solution u 6 p + [a, (3] of = u. 

While this is reported to provide tight numerical enclosures at reasonable cost 
for computations at machine precision even in the case of nonlinear equations in 
many variables, we are not aware of any asymptotic tightness result of the kind 
we are interested in in this paper. In fact, it is not entirely clear to us under 
which conditions methods of this kind are guaranteed to produce arbitrarily tight 



enclosures. (Note however that van der Hoeven (20071 states initial results in this 



direction.) Neither do we know how to use them to bound tails of D-finite functions 
on their whole disk of convergence. 

And yet, D-finite functions may be evaluated to an absolute precision 10~™ in 
softly linear time n(logn)°^ by computing truncations of their Taylor series by bi- 
nary splitting. Numerical analytic continuation based on this technique then allows 



to obtain values of these functions at any point of their Riemann surfaces (Chud 



novsky and Chudnovsky 1988 §5). Applications include the numerical computa- 



tion of monodromy matrices of linear differential equations with polynomial coef- 
ficients. In this context, one benefit of the language of majorant series is that a 
single majorant encodes both bounds on the values and truncation orders for all 
elements of a basis of the local solutions of the differential equations as well as 
their derivatives — all of which are useful to control errors in the numerical analytic 
continuation process. 

Excluding degenerated cases, the number of terms of the series to take into 
account is Xn+o(n), where A depends on the location of the evaluation point relative 
to the singularities of the function, or Oinj log n) in the case of entire functions. 
The tightness result of Theorem [l] translates into the fact that the number N of 
terms that get computed is indeed of that order, while most existing methods for 
computing bounds of tails of D-finite series seem to ensure only N = O(n). This 
in turn improves the complexity of the algorithm by a constant factor. 

The subpackage of gfun mentioned above contains high-precision numerical eval- 
uation and analytic continuation routines based on this strategy. They rely on the 
code computing bounds for accuracy control. These numerical evaluation facilities 
are exported to the DDMF. 
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5.3. Experiments. In Table[T] we report on experiments concerning the tightness 
of the bounds for truncating Taylor series expansions of a few common elementary 
and special functions. Each column label actually stands for a differential equation 
that annihilates the given function (with suitable initial values), and an evaluation 
point smaller in absolute value than the dominant singularity of the differential 
equation. Each internal cell shows the truncation order computed by the algorithm 
from this data for a specific accuracy requirement, and compares it to the minimal 
correct answer, computed by exhaustive search. For instance, the column "erf(l) 2 " 
corresponds to the evaluation at z = 1 of the function u(z) — erf(z) 2 represented 
as the unique solution of 

(2 + 8z 2 )u'(z) + 6zu"(z) + u"'(z), u(Q) = 0, u'(0) = 0, u"(0) = 8 . 

7T 

Using a majorant series for u, our algorithm determined that |u ; i9o(l) — tt(l)| < 
10 00 , but it happens that only the first 163 of these 190 terms are really necessary. 
It can be seen that the bounds we compute do not stray too far from the optimal 
values. 

We consider three cases, corresponding to the three main types of asymptotic 
behaviours that the coefficient sequence of a convergent D-finite series may exhibit, 
characterized (in generic cases) by the nature of the dominant singularities of the 
differential equation: regular singularities (k = = T with the notations of the 
previous sections), irregular singularities at finite distance (k = 0,T > 0), or at 
infinity (k < 0). (Irregular singularities with n > correspond to divergent power 
series, and a differential equation whose only singularity is a regular singular point 
at infinity has only polynomial solutions. The examples of the second set all involve 
right composition by rational functions because it is unusual to study differential 
equations with more than two irregular singular points, and those are usually taken 
to be oo and 0.) 

For each of these, the last three columns illustrate how the truncation orders and 
the bounds vary as \z\ approaches the radius of convergence of the series. Note that 
high-order Taylor expansions at are not the best way to compute numerical values 
of D-finite functions for such z: the growth of the truncation orders (both optimal 
and computed) can be got around by using several steps of analytic continuation 
along a broken- line path from to z ( Chudnovsky and Chudnovsky| 1987| §4). 



The example of Si(z) has an interesting feature: the origin is a regular singular 
point of the differential equation mentioned in Example |TJjd|) , but Si(z) may never- 
theless be defined by simple initial values at origin, so that our algorithm applies 
without any adjustment. 

Finally, here is a nontrivial "non-generic" example where our method fails to 
produce a tight bound. 



Example. In his proof or the irrationality of £(3), Apery (19791 introduces two 
sequences (a„) and (&„) such that u n = b n — C(3)a„ satisfies the (minimal-order) 
linear recurrence relation 

(n+2) 3 u n+2 = (2n+3)(17n 2 +51n+39) u n+1 -(n+l) z u n , u = -£(3), u x = 6-5C(3). 
Applied to this recurrence relation, Algorithm [5] determines that 

| tin | < 1-21 (n 2 + 3n + 2) (17 + 12^2)" (where (17 + 12\/2) ~ 33.97) 
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This bound is asymptotically tight for both a n and b n , but the whole point of 
Apery's proof is that b n — ((3)a n — > fast as n — > oo. 

Acknowledgements. We thank Moulay Barkatou, Nicolas Brisebarre, Sylvain Chevil- 
lard and Nicolas Le Roux for interesting discussions or comments on earlier versions 
of this work that have led to improvements, and an anonymous referee for spotting 
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Regular dominant singularity 
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Finite irregular dominant singularity 
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Dominant singularity at infinity 
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Table 1. Computed/minimal required number of terms of the 
Taylor expansion of a D-finite function to approximate this func- 
tion to a given absolute precision. In this table, ip is the solution 
of the spheroidal wave equation (1 — z 2 ) ip"(z) — 2(6 — l)zip'(z) + 
(c — 4qz 2 )tp(z) = given by the choice of parameters and initial 
values b = 1/2, q = 1/3, c = 1, ip(0) = 1, ip'(0) = 0; Ai and Bi 
denote the Airy functions; erf stands for the error function and Si 
for the integral sine. 



